module mod_matrix
  implicit none
  
  !
  ! lapack-2D matrix
  !
  !
  ! (c) Nicolaas E. Groeneboom
  !     leuat@irio.co.uk
  !     http://www.irio.co.uk
  !

 type, public :: t_matrix
     real*8, dimension(:,:), allocatable :: val
     integer :: size
     real :: determinant
  end type t_matrix

 
contains  

subroutine matrix_initialize(matrix, size)
    implicit none
    type(t_matrix), intent(inout) :: matrix

    integer, intent(in) :: size

    if (allocated(matrix%val) .and. (matrix%size .ne. size)) deallocate(matrix%val)
    matrix%size = size
    if (.not. allocated(matrix%val))  allocate(matrix%val(size,size))
    matrix%val = 0.0
    

end subroutine matrix_initialize

subroutine matrix_print(matrix)
  type(t_matrix), intent(inout) :: matrix
  integer ::i,j
  do i = 1, matrix%size
     write(*,*) matrix%val(i,:)
  end do
  

end subroutine matrix_print

subroutine matrix_invert(from_matrix, to_matrix) 
  external DGETRF
  external DGETRI
  
  type(t_matrix), intent(inout) :: from_matrix
  type(t_matrix), intent(inout) :: to_matrix
  integer :: info, n, i
  real, dimension(:), allocatable :: work  ! work array for LAPACK
  integer, dimension(:) , allocatable :: ipiv   ! pivot indices
  
  n = from_matrix%size

  allocate(work(n*n))
  allocate(ipiv(n))
  ipiv = 0.0
  work = 0
 
  to_matrix%val = from_matrix%val
  
 
  call DGETRF(n, n, to_matrix%val, n, ipiv, info)
  
  !if (info /= 0) then
     !stop "Matrix is numerically singular!"
  !end if
  if (info==0) then
  to_matrix%determinant = 1
  do i=1, n
     to_matrix%determinant = to_matrix%determinant*to_matrix%val(i,i)
  end do
  !write(*,*) "Determinant ", to_matrix%determinant

  call DGETRI(n, to_matrix%val, n, ipiv, work, n, info)

  if (info /= 0) then
     stop 'Matrix inversion failed!'
  end if

  end if

  deallocate(work)
  deallocate(ipiv)

  
end subroutine matrix_invert

 
subroutine matrix_deallocate(matrix)
    implicit none
    type(t_matrix), intent(inout) :: matrix
    if (allocated(matrix%val)) deallocate(matrix%val)
    matrix%size = -1
  end subroutine matrix_deallocate



function matrix_trace(matrix) 
    implicit none
    type(t_matrix), intent(inout) :: matrix
    integer :: i
    real*8 :: matrix_trace

    matrix_trace = 0
    do i =1, matrix%size
       matrix_trace = matrix_trace + matrix%val(i,i)
    end do
    
  end function matrix_trace
 
subroutine matrix_identity(matrix)
    implicit none
    type(t_matrix), intent(inout) :: matrix
    integer :: i
    matrix%val = 0
    do i =1, matrix%size
       matrix%val(i,i) = 1.0
    end do
  end subroutine matrix_identity


SUBROUTINE gaussj(a,n,b) 
INTEGER n,NMAX 
REAL*8 a(n,n),b(n) 
PARAMETER (NMAX=50) 
INTEGER i,icol,irow,j,k,l,ll,indxc(NMAX),& 
        indxr(NMAX),ipiv(NMAX) 
REAL*8 big,dum,pivinv 
do j=1,n 
  ipiv(j)=0 
end do 
do i=1,n 
  big=0. 
  do j=1,n 
    if(ipiv(j)/=1) then 
      do k=1,n 
        if (ipiv(k)==0) then 
          if (abs(a(j,k))>=big) then 
            big=abs(a(j,k)) 
            irow=j 
            icol=k 
          endif 
        else if (ipiv(k)>1) then 
          pause 'singular matrix in gaussj' 
        endif 
      end do 
    endif 
  end do 
  ipiv(icol)=ipiv(icol)+1 
  if (irow/=icol) then 
    do l=1,n 
      dum=a(irow,l) 
      a(irow,l)=a(icol,l) 
      a(icol,l)=dum 
    end do 
    dum=b(irow) 
	b(irow)=b(icol) 
	b(icol)=dum 
  endif 
  indxr(i)=irow 
  indxc(i)=icol 
  if (a(icol,icol)==0.) pause 'singular matrix in gaussj' 
  pivinv=1./a(icol,icol) 
  a(icol,icol)=1. 
  do l=1,n 
    a(icol,l)=a(icol,l)*pivinv 
  end do 
  b(icol)=b(icol)*pivinv 
  do ll=1,n 
    if(ll/=icol) then 
      dum=a(ll,icol) 
      a(ll,icol)=0. 
      do l=1,n 
        a(ll,l)=a(ll,l)-a(icol,l)*dum 
      end do 
      b(ll)=b(ll)-b(icol)*dum 
    endif 
  end do 
end do 
do l=n,1,-1 
  if(indxr(l)/=indxc(l)) then 
    do k=1,n 
      dum=a(k,indxr(l)) 
      a(k,indxr(l))=a(k,indxc(l)) 
      a(k,indxc(l))=dum 
    end do 
  endif 
end do 
END SUBROUTINE gaussj 


end module mod_matrix


